erweiterte Näherung zu Pythagoras

← Zurück zum Forum

Dr. Irata 02.05.2025 14:38

... der Pythagoras läßt mich einfach nicht los!

In meinen langen Erläuterungen und Versuchen einen schönen Raycaster auf dem A8 zu bauen (ich denke das Projekt werde ich wieder im Herbst aufnehmen), bin ich zwingend abhängig auf die Berechnungen von Entfernungen. In Assembler ist das gar nicht so einfach, denn quadrieren und Wurzel ziehen bedarf besonderer Routinen (die ich ja später hier alle hergeleitet habe). Diese Routinen existierten damals so nicht und ich war auf eine ganz gute Näherung angewiesen, die ich dann auch im Netz finden konnte - in den 50 igern und 60 igern gerade bei den einfachen Rechnern waren diese Näherungen notwendig und einige Mathematiker echt erfinderisch. Die Näherung für c^2 = a^2+b^2 also c = sqrt (a^2+b^2) ist ja: c = a+b/2 für a>b Die Näherung ist echt gut und super schnell. Leider ist die an bestimmten Stellen fehlerhaft. Diesen Fehler kann man aber ganz gut ausgleichen: Tendenziell ist die Näherung immer zu groß. Wenn man jetzt aber vom Ergebnis c die erste Stelle von a+1 subtrahiert solange b>10 ist, dann ist das Ergebnis wirklich nahezu perfekt. Beispiel: a=40 , b=30 --> c=40-5+30/2 = 50 (mathematisch 50) Beispiel: a=35 , b=25 --> c=35-4+25/2 = 43 (mathematisch 43) Beispiel: a=72, b=62 --> c=72-8+62/2 = 95 (mathematisch 95)

Dr. Irata 02.05.2025 17:48

Also: c = sqrt (a^2+b^2) wenn a>=b wenn b>=10 dann c=a-(a/10)+b/2 wenn b<10 dann c=a+b/2

wenn b>a wenn a>=10 dann c=b-(b/10)+a/2 wenn a<10 dann c=b+a/2

und jetzt eine Aufgabe in die "Gemeinde":

Wer kann das in möglichst schnellen Assembler-Code umsetzten? Bin gespannt auf eure innovativen Lösungen...

RhoSigma 02.05.2025 18:39

Äh, hab zwar keinen Code für Dich, ABER warum Basis 10? Die division durch 10 ist lahm. Nimm doch 16. Oder 8?

Dr. Irata 02.05.2025 18:51

... alle guten Ideen sind hier willkommen. Ich brauche ja die erste Ziffer der Zahl a - also geteilt durch 10. Das ist blöd in Assembler. Durch 8 teilen passt bestimmt auch noch für die Näherung - oder im Dezimalmodus den vorderen Nibble extrahieren... Es muss halt schnell sein! Wenn ich durch 8 teile, dann fällt sogar das +1 weg ... also ist noch besser. Das funktioniert bis 70 sehr gut. Viel größere Entfernungen werde wir am A8 bei den möglichen Spielen nicht brauchen. Sehr gut!!

Kveldulfur 02.05.2025 19:07

Hallo!

        ORG $2000

DIV10TAB
        .BYTE 0,0,0,0,0,0,0,0,0,0
        .BYTE 1,1,1,1,1,1,1,1,1,1
        .BYTE 2,2,2,2,2,2,2,2,2,2
        .BYTE 3,3,3,3,3,3,3,3,3,3
        .BYTE 4,4,4,4,4,4,4,4,4,4
        .BYTE 5,5,5,5,5,5,5,5,5,5
        .BYTE 6,6,6,6,6,6,6,6,6,6
        .BYTE 7,7,7,7,7,7,7,7,7,7
        .BYTE 8,8,8,8,8,8,8,8,8,8
        .BYTE 9,9,9,9,9,9,9,9,9,9
        .BYTE 10,10,10,10,10,10,10,10,10,10
        .BYTE 11,11,11,11,11,11,11,11,11,11
        .BYTE 12,12,12,12,12,12,12,12,12,12
        .BYTE 13,13,13,13,13,13,13,13,13,13
        .BYTE 14,14,14,14,14,14,14,14,14,14
        .BYTE 15,15,15,15,15,15,15,15,15,15
        .BYTE 16,16,16,16,16,16,16,16,16,16
        .BYTE 17,17,17,17,17,17,17,17,17,17
        .BYTE 18,18,18,18,18,18,18,18,18,18
        .BYTE 19,19,19,19,19,19,19,19,19,19
        .BYTE 20,20,20,20,20,20,20,20,20,20
        .BYTE 21,21,21,21,21,21,21,21,21,21
        .BYTE 22,22,22,22,22,22,22,22,22,22
        .BYTE 23,23,23,23,23,23,23,23,23,23
        .BYTE 24,24,24,24,24,24,24,24,24,24
        .BYTE 25,25,25,25,25,25     


VarA        .BYTE 50
VarB        .BYTE 30
VarC        .BYTE 00

Calc:       LDA VarA        ; A>=B?
        CMP VarB
        BCS AgtB
        JMP BgtA

AgtB:       LDA VarB        ; B>=10?
        CMP #10
        BCS Bgt10

        LDA VarB
        LSR         ; b/2
        CLC
        ADC VarA
        STA VarC

        RTS

Bgt10:      SEC
        LDX VarA    
        TXA
        SBC DIV10TAB,X      ; a/10
        STA VarC        ; VarC = a-(a/10)

        LDA VarB        ; b/2
        LSR

        CLC
        ADC VarC        ; a-(a/10)
        STA VarC        ; VarC = b/2 + a-(a/10) 

        RTS

BgtA:       LDA VarA        ; a>=10?
        CMP #10
        BCS Agt10

        LDA VarA
        LSR         ; a/2
        CLC
        ADC VarB
        STA VarC

        RTS

Agt10:      SEC
        LDX VarB    
        TXA
        SBC DIV10TAB,X      ; b/10
        STA VarC        ; VarC = b-(b/10)

        LDA VarA        ; a/2
        LSR

        CLC
        ADC VarC        ; b-(b/10)
        STA VarC        ; VarC = a/2 + b-(b/10) 

        RTS

Start:      JSR Calc
        BRK

        RUN Start

Wenn nur 8 Bit benötigt werden, kann man geteilt 10 am Besten in einer Tabelle realisieren. Ich hoffe der Code hat keine Flüchtigkeitsfehler 🙂

Grüße Janko

Meine Projekte findest Du hier...

Dr. Irata 02.05.2025 19:50

Hi Janko, sehr super und danke. Ich habe das aber nochmal geprüft. Nimm mal nicht durch 10 sondern durch 8. Dann kann man das +1 für den zu subtrahierenden Wert weglassen. Das +1 hatte ich ganz oben erwähnt, in der Beschreibung für den Code aber vergessen. Wenn man aber statt durch 10 nur durch 8 teilt, dann fällt praktischerweise das +1 weg.

Wieviele Zyklen braucht das dann?

Dr. Irata 02.05.2025 20:30

Janko hat mir gerade geschrieben, daß das so etwa 34-51 Zyklen braucht. Das ist echt super schnell - und läßt sich sicher noch ein wenig optimieren. Wer schafft es schneller?

Diese Näherung ist wirklich gut, schwächelt nur etwas bei b<10 ... da werde ich mal schauen, ob man das noch für den Bereich besser schafft.

RhoSigma 02.05.2025 23:01 
        ORG $2000

VarA        .BYTE 50
VarB        .BYTE 30
VarC        .BYTE 00

Calc:
        LDA VarA        ; A>=B?
        CMP VarB
        BCS AgtB

BgtA:
        CMP #10     ; a>=10?
        BCS Agt10

        LSR         ; a/2
        CLC
        ADC VarB
        STA VarC

        RTS

AgtB:
        LDA VarB        ; B>=10?
        CMP #10
        BCS Bgt10

        LSR         ; b/2
        CLC
        ADC VarA
        STA VarC

        RTS

Bgt10:
        LDA VarA    
        LSR
        LSR
        LSR         ; a/8
        SEC     
        SBC VarA        ; (a/8)-a
        EOR 0xFF
        TAX
        INX         ; a-(a/8)
        STX VarC        ; VarC = a-(a/8)

        LDA VarB        ; b/2
        LSR

        CLC
        ADC VarC        ; a-(a/8)
        STA VarC        ; VarC = b/2 + a-(a/8) 

        RTS


Agt10:
        SEC
        LDX VarB    
        TXA
        SBC DIV10TAB,X      ; b/10
        STA VarC        ; VarC = b-(b/10)

        LDA VarA        ; a/2
        LSR

        CLC
        ADC VarC        ; b-(b/10)
        STA VarC        ; VarC = a/2 + b-(b/10) 

        RTS

Start:
        JSR Calc
        BRK

        RUN Start

So, Variante mit /8 und 2 Instruktionen gespart, ungetestet natürlich ... Ohne Tabelle ist Div/8 sogar langsamer, aber man könnte auch die Tabelle benutzen und die Zyklen wieder einsparen. Dafür ist es halt kürzer.

Dr. Irata 03.05.2025 08:57

ChatGPT berechnet daraus zwischen 26 und 60 Taktzyklen - je nach Werte und Bedingungen. Das ist auf jeden Fall schön schnell!

Dr. Irata 03.05.2025 09:10

Bei Jankos Ansatz müsste man ja nur die Tabellenwerte um den Wert 1 erhöhen. Die Taktzyklen bei Janko sind zwischen 31 und 47 - das sind im Mittel aller 4 Fälle 39,75 Taktzyklen. Der Ansatz von RhoSigma liegt bei 41,5 Taktzyklen - also im Mittel etwas langsamer.

Ich werde hier bei Zeiten noch einen etwas anderen Ansatz mit direkten Look-Up-Tabellen posten, der etwas genauer ist und noch etwas schneller sein könnte (ich hoffe auf stabil 31 Zyklen).

RhoSigma 03.05.2025 21:04

ChatGPT-4o schlägt noch diese Näherung vor:

c≈a+b− (a⊕b)/4

Die wäre sehr effizient zu implementieren:

    org $2000

        ; Zero-page variables for extra speed
a   .byte 5
b   .byte 3
t   .byte 0

LENGTH:
        LDA a
        EOR b
        LSR
        LSR
        ADC #0 ; round 0.5 up
        STA t  ; t = (a XOR b) / 4

        LDA a
        CLC
        ADC b ; c = a + b

        SEC
        SBC t
        RTS     ;  final result in A

MAIN:
    JSR LENGTH
        BRK       ; End of program

Die Genauigkeit ist aber nicht besonders gut. Geringfügig besser wird es, wenn du /4 rundest, also

c≈a+b− ((a⊕b) + 2)/4

Macht im Code ADC #0 mehr (nach dem 2. LSR) ​

Dr. Irata 03.05.2025 21:41

die ist allerdings deutlich länger als c = a + b/2 , a>b

LarsImNetz 04.05.2025 18:37

Ich nutze in meinem aktuellen Spiel die einfache Variante a+b/2 für a<b Genauer:

value = y<x ? x + y/2 : y + x/2

y und x sind vom Typ INT8 also [-128;127] ich nutze aber nur [-64;63] das kostet

Das reicht völlig aus, jedenfalls bei mir.

Als Tipp: Probiere es doch erstmal mit dieser Version aus, ist die Genauigkeit nicht ausreichend, dann nimmst halt etwas besseres, dann wird es ein wenig langsamer, na und. Du rechnest ja nicht nur Hypotenusen aus, sondern das Programm macht noch jede Menge mehr.

Wenn es richtig auf Speed ankommt, kannst Du für die paar Berechnungen den Code vielleicht Inline einfügen, das spart pro Aufruf 12 Zyklen!

JM2C.

Dr. Irata 04.05.2025 19:50

du meinst aber sicherlich c = a + b/2 für a>b

RhoSigma 04.05.2025 20:42

Hier ist noch eine Variante, die im Punkt Fehlerabweichung recht ähnlich ist wie Deine ursprüngliche Variante (mit der Korrektur für min >= 10).

c = max(a,b) + (min(a,b)+1)/2 - (min(a,b)+4)/8

length:
    LDA b      ; Load b into A                   ; 2 cycles
    CMP a      ; Compare b with a                ; 3 cycles
    BCS a_is_max ; Branch if a >= b (C is set if a >= b) ; 2 cycles (not taken), 3 cycles (taken)

b_is_max:
; Calculate b + a / 2 - a / 8
    LDA a      ; Load a into A                   ; 3 cycles
    LSR        ; A = a / 2                       ; 2 cycles
    LSR        ; A = a / 4                       ; 2 cycles
    LSR        ; A = a / 8                       ; 2 cycles
    ADC #0     ; (a + 4)/8                       ; 2 cycles
    STA temp   ; store temp                      ; 3 cycles

    LDA a      ; Load a into A                   ; 3 cycles
    LSR        ; A = A / 2                       ; 2 cycles
    ADC b      ; A = b + ((a+1) / 2)             ; 3 cycles
    SEC                                          ; 2 cycles
    SBC temp   ; final result in A               ; 3 cycles
    RTS

a_is_max:
; Calculate a + b / 2 - b / 8
    LSR        ; A = b / 2                       ; 2 cycles
    LSR        ; A = b / 4                       ; 2 cycles
    LSR        ; A = b / 8                       ; 2 cycles
    ADC #0     ; (b + 4)/8                       ; 2 cycles
    STA temp   ; store temp                      ; 3 cycles

    LDA b      ; Load b into A                   ; 3 cycles
    LSR        ; A = A / 2                       ; 2 cycles
    ADC a      ; A = a + ((b+1) / 2)             ; 3 cycles
    SEC                                          ; 2 cycles
    SBC temp   ; final result in A               ; 3 cycles
    RTS

34-36 Zyklen, wenn ich mich nicht verrechnet habe und a, b, temp ZP sind. Man kann noch 2 Zyklen sparen, wenn man auf die SEC vor den SBC temp verzichtet, das macht keinen grossen Unterschied (de facto zieht es IMMER eine 1 im interessanten Bereich ab).

Ausserdem produziert diese Variante ein "stetiges" Resultat, Deine hat Ausreisser, siehe Bilder.

Figure_my.png Figure_orig.png

Anhänge:
Dr. Irata 04.05.2025 21:05

... echt super ... ich hatte diese blöden "Ausreißer" schon wahrgenommen aber noch nicht so richtig in den Griff bekommen! Echt gut!!

FlorianD 04.05.2025 21:22

ich kann zwar nichts beitragen, weil mein Mathe nach "+ - * /" endet, aber bin völlig fasziniert, wie ihr hier ein Problem seziert und löst!! Weiter so!

RhoSigma 04.05.2025 22:10

Dr. Irata schrieb: ... echt super ... ich hatte diese blöden "Ausreißer" schon wahrgenommen aber noch nicht so richtig in den Griff bekommen! Echt gut!!

Leider geht das Lob an ChatGPT, dass mir den python Rahmen "programmiert" hat. War echt das erste Mal, das da was sinnvolles rauskommt 🙂 hab ihm auch Dein Problem versucht schmackhaft zu machen, aber der phantasiert immer irgendwelche nichtexistierenden 6502 Befehle zusammen und ich kriege es nicht ausgetrieben: er sagt immer "you are right, ADC X is a non-existent instruction" und ersetzt das dann mit ADC Y 😀 Offenbar gab es in Stackoverflow mehr python code wie 6502 Assembler "zu erlernen"... Wenn gewünscht, poste ich hier den Python Code für weitere Experimente...

Dr. Irata 05.05.2025 07:27

... tatsächlich habe ich mich nach Anregung von Janko an diesem WE zum ersten Mal mit ChatGPT beschäftigt und mich direkt bei der deutschen Variante kostenfrei angemeldet. Hier bin ich etliche Dinge mit Mike (ich habe ihn Mike genannt) durchgegangen und war echt erstaunt, wie gut das funktioniert. Mike hat auch wirklich gute Ideen zu allen Fragen der Assemblerprogrammierung unter 6502. Mike konnte mir auch weiterhelfen bei der Reparatur meines Boseboards... insgesamt ist das natürlich eine lernfähige PseudoKI - sie sucht sehr zielgerecht und lernend das gesamte Internet noch Lösungen auf deine Frage durch und baut dann eine halbwegs sinnvolle Antwort zusammen. Wenn in ein paar Jahren die Quantencomputer serienreif werden, dann bin ich mir ganz sicher: Viele Dinge werden nachher nicht mehr so sein, wie vorher!!!